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vQ [ Nanostructures, like periodic arrays of scatters or low-index gratings, are used to improve the light outcoupling 
^-H 1 from organic light-emitting diodes (OLED). In order to optimize geometrical and material properties of such 
structures, simulations of the outcoupling process are very helpful. The finite element method is best suited 
' ly^' , for an accurate discretization of the geometry and the singular-like field profile within the structured layer and 
O ' the emitting layer. However, a finite element simulation of the overall OLED stack is often beyond available 
-j-J I computer resources. The main focus of this paper is the simulation of a single dipole source embedded into a 
p~^' twofold infinitely periodic OLED structure. To overcome the numerical burden we apply the Floquet transform, 
so that the computational domain reduces to the unit cell. The relevant outcoupling data are then gained by 
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ABSTRACT 



Yi ' inverse Flouqet transforming. This step requires a careful numerical treatment as reported in this paper. 
c/3 Keywords: organic light emitting diodes, light extraction. Green's tensor, Floquet transform 

^ ■ 
P^' 1. INTRODUCTION 

Figure [1] shows a simplified OLED structure. It essentially consists of a layered medium stack. Light is generated 
in a slim emitting layer within the organic semiconductor material. Often, the metallic cathode also serves as 
an optical back reflector. The goal is to extract as much of the generated light as possible into the superstrate. 
Full extraction is inhibited due to the presence of lossy materials and by the trapping of light into waveguide 
Z^ ' modes by total reflection. The waveguide modes travel in horizontal direction and are therefore lost for emission. 
^^ , Scattering particles are commonly used to disrupt the propagation of the waveguide modes; the light, gradually 
_il ' scattered by the particles, can then leave the OLED device. Unfortunately, periodically arranged scatterers are 
(•~^ , in general not able to guarantee full light extraction even not for transparent materials: Light is still trapped in 
C^ ■ Bloch modes (except for frequencies within the band gap of the twofold photonic crystals). A proper design of 
^"^ , the periodic arrangement is therefore of major importance for an efficient OLED. 

^ , The simulation of light extraction properties is a numerical challenge for the following reasons: 



1. A wavelength scan over the entire visible spectrum is needed. 

2. Many light emitters have to be simulated at different positions within the structured OLED. 

3. Metals gives rise to the presence of plasmons with singular field profiles near metallic edges or corners. 

4. Realistic material data are only given experimentally. A numerical dispersion model can be costly to 
implement. 



Further author information: (Send correspondence to Lin Zschiedrich, E-mail: lin.zschiedrich@jcmwave.com .) 





X 

emitters 
scatterers 




^^ r \ * ' y superstrate 

anode 
organic 

organic 

cathode 
substrate 

Figure 1: Simplified OLED structure. Light is generated in a shm emitter layer and radiated into the upper 
half-space. Tiny scatterers are used to increase the light extraction efficiency. 

5. The computational domain must be sufficiently large to suppress numerical truncation errors. 

In this paper we focus on the finite element method (FEM) in the frequency domain. This method allows for 
an efficient and accurate discretization of the geometry as well as an automatic mesh adaption for an accurate 
resolution of the highly nonuniform field profiles. Since we apply FEM in frequency domain each wavelength is 
treated separately. However, multiple dipole source positions can be computed efficiently in one sweep. This is 
because the finite element method is chiefly limited by the direct sparse matrix solver, which can be re-used for 
different source terms. 

In this paper we discuss a method which allows to simulate an isolated source embedded in a twofold periodic 
arrangement without the need to use a large computational domain. By means of the Floquet transform, see 
Kuchment in Gao et al. [U pp. 207], the original problem posed on the entire periodic space is mapped to 
a bundle of Bloch-periodic problems posed on the unit cell of the periodic structure. These Bloch-periodic 
problems can be solved with a tremendous reduction of memory requirements. The price we have to pay is the 
inverse Floquet transform which is an integration over the Brillouin zone of the reciprocal lattice space. Using 
adaptive integration techniques together with a straightforward parallelization we show that this can be done 
with reasonable numerical costs. The idea to numerically employ the inverse Floquet transform goes back to 
Wilcox, Botten, McPhedran et al.^ There, the motivation was the computation of defect modes in photonic 
crystals. It has been shown that the inverse Floquet transform is still numerically feasible even close to the 
band-edge of the photonic crystal. 

The paper is organized as follows: In Section[5]we settle the basic concepts for modelling light extraction from 
a light emitting diode. Then we explain how the finite element method accurately deals with singular sources 
such as point dipoles (Section|3]). Section|3]introduces the Floquet-transform techniques. The final section covers 
numerical concepts and examples. 



2. MAXWELL'S EQUATIONS: LIGHT EMISSION MODEL 

We consider Maxwell's equations in the frequency domain. That is, we assume a time-harmonic dependency of 
the electromagnetic field, i.e.. 
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E(r,i) = Re^EC 

H(r,i) = Re('H(r,w)e~*'^*V 
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Figure 2: Computed efficiencies for a planar OLED stack. Left: horizontally oriented dipole. Right: vertically 
oriented dipole. The lower substrate consists of an infinite silver cathode (Drude model) followed by layers 
(n: refractive index (i:thickness) ni = 1.75, di = 80nm; 712 = 1.75,^2 = lOOnm; n^ = 1.8,^3 = 200nm; 
n4 = 1.9, ^3 = 600nm, and an glass superstrate with n = 1.5. The dipoles are placed between the first and 
second layer. This non-trivial example demonstrates the fidelity of the FEM approach. 



and accordingly for the source current J. The constitutive relations are given in the form 

i>{r,uj) = e(r,a;)E(r,w), 
B(r,w) = /i(r,u;)H(r,w), 
J(r,w) = cr(r,w)E(r,w) + J.,(r,a;), 



(la) 
(lb) 
(Ic) 



with the material tensors permittivity, e, permeablity, /i, and conductivity, a. 3i{r,uj) is the impressed current 
density. For simplicity we will henceforth drop the hats, will use J for the the impressed source, and we introduce 
the complex permittivity tensor e — e + ia juj. Then, Maxwell's equations can be cast into a second order form 
for the electric field. 



V X ^"^V X E(r) - u?£¥L{v) ^ iwJ(r). 



(2) 



In an OLED simulation, the major quantity of interest is the extraction efficiency. This is defined as the quotient 
of the power Prad, radiated into the superstrate, and the total emitted power Ptot of the source. 
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More precisely, the radiated power is the integrated power fiux through an infinitely far upper hemisphere, 
that is. 
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where Sr^^ is the upper hemisphere with radius R and n is normal vector. The last equality holds true, because 
the far field satisfies the Silver-Miiller radiation condition, see Monk [3, p. 226]. 

To compute the total emitted power Ptot we regard a domain f2 containing the source. Then, the emitted 
power is the energy lost in this domain plus the net power fiux across the boundary f2. 



^tot = / ^Re (E X H) d5 + / ^E • aEdV. 



This expression can be simplified to a linear functional for the electric field. To see this, we multiply Maxwell's 
equations ^ with E and integrate over J7: 



in 
Partially integrating yields 



/ E- V X /i^^V X E-tj^E-eEdV"^ iuj J^^E ■ J dV. 
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V X E-z^^-'V X E-cj^E-eEdV 
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I (E X ^"^V X E) -n dS* = iw/j-jE-JdK 
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Recalling that Im (e) = a/u and taking the imaginary part on both sides gives 



Ptot = -^ /Re(E-J) dV. 
^ Jn 



3. MODELLING DIPOLE SOURCES WITH FINITE ELEMENTS 

For a dipole source at position r' the impressed electric current is modeled as a delta distribution, J(r) = p(5(r— r') 
with given dipole moment p. The regularity of E is poor and a direct finite element discretization of E suffers from 
a slow convergence. To cure this we use the subtraction approach, see Awada et al.,SjWolterfP and Zschicdrich.S) 

The idea behind the subtraction approach is to determine an analytically available singular field Es which 
already comprises the singular part of E at the dipole position. A natural candidate is the homogeneous Green's 
function E^, that is 

V X ^-^V X E,(r) - uj'^ea'E.siv) = iu^5{v - r'), 

with a constant material background ed = e(r') and /i^ = /^(r') as given at the dipole position. Now, we split 
the field E into the singular field E^ and a correction field Ec that is E = Es + Ec. Inserting into Maxwell's 
equations ([2]) yields 

V X /i-iy X (E, + Ee) (r) - Lo^e (E, + E,) (r) = 
V X /i^^V X Ec(r) - u'^e'Ecir) + V x ^"^V x Es(r) - uj'^e'Esir) = 
V X /i-iy X Ee(r) - w^eEeW + V x (^"^ - ^-i)V x Es(r) - ^^(e - ed)E,(r) + 

V X Aid ^V X Es(r) - Lo'^Ed + Es(r) =iwp(5(r - r') . 

— zcjp5(r— r') 

Hence, as desired, the singular source terms on both sides cancel out. When rearranging the terms we end up 
with Maxwell's equations for the correction field Ec only, 

V X ^"^V X Ec(r) - w^eEc(r) = 

- V X (//-i - /i-i)V X Es + w2(e - erf)Es. 

The analytically given right hand side is equal to zero in a vicinity of the dipole position. Hence, this equation 
for the correction field Ec is well suited for an accurate discretization with finite elements. (Even a jump in the 
permeability is allowed, as in the variational form the most left V x — operator on the right hand side can be 
applied on the test function by partial integration.) 
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Figure 3: An OLED stack with an periodically structured cathode. Only finite layers are plotted. The material 
stack is given in the table. The dipoles are placed between the second and third layer (red plane) 

Remark 

For simplicity we have assumed so far that the material background is homogeneous in a vicinity of the dipole. 
The case of a dipole placed on or near a material interface can be treated in the same way: As the singular field 
Es, one the chooses the solution of the dipole source for a two layer material, which is available quasi-analytically, 
see Paulus.'^ 

It remains to evaluate the radiation efficiency ?7iad in an accurate manner for the dipole case. The total 
radiated power can be computed as usual from the far field data. But the expression ^ for the total emitted 
power calls for a delicate mathematical justification as it involves the integral of the delta distribution with a 
singular field. Fortunately, one can show that this expression is properly defined when the dipole is placed in a 
lossless background. Then, the total emitted power is given by 



Ptot = -^Re(E(r')-p) 



(4) 



Again, we evaluate the expression by using the splitting E = E^ + Ec. Re (Es(r') • p) is analytically available, 
whereas Ec is smooth and can be evaluated within the finite element framework with high precision. 

To validate the FEM method we compare the results for a planar OLED stack with the quasi-analytic solution 
as obtained by Fourier expansion techniques, c.f., Paulus et al.^ This is a non-trivial test for the FEM approach, 
since only the homogeneous dipole solution was used as the singular field. Figure [5] shows a great agreement of 
the 3D numerical solution (FEM) with the exact (analytic) result. 

4. PERIODIC GEOMETRIES 

We now address the case of a twofold periodically structured device, that is 

e(r + ai/2) =e(i-), 
/i(r + ai/2) =/i(r), 

with grid vectors ai, a2 in the xy— plane. 

In the following, let 1 = [^i; I2] G N^ denote an integer vector and we use the notation a = [ai, a2] (in matlab 
style). We call a source field Bloch-periodic, when it satisfies 



J(r + a-l) = e*Ba-ij(i.) 

with a Bloch-phase vector ke G R^. The corresponding electric field is Bloch-periodic as well. This allows to 
restrict the computation onto an unit cell by imposing Bloch-periodic boundary conditions. 



However, the main focus of this paper is the simulation of a single dipole source embedded into a periodic 
arrangement. Since this single dipole source is not Bloch-periodic, the entire space is in principle needed as 
computational domain. However, in the sequel we will explain that the computational domain can be reduced to 
the unit cell by using the Floquet transform. For any sufhciently decaying source term we perform the Floquet 
transform on J : 



^e*S'^'j(r-a-l). 



IGF 



One readily verifies that Jkg is Bloch-periodic with phase vector ke- Hence the corresponding Bloch-periodic 
solution Eke can be computed on a unit cell. The solution E for the orginal source term J is then obtained by 
the inverse Floquet transform. For doing this we introduce the reciprocal lattice vectors bi, b2 satisfying 



[bi,b2]'^ • [ai,a2] = 27r 



1 
1 



with bi, b2 perpendicular to the plane spanned by ai, a2. Again we write b — [bi, b2] and define the Brillouin 
zone 

BZ = {kB =b-r Ire [0,1] x [0,1]}. 

Integrating Jk^ over the Brillouin zone reproduces the original source field J : 



I Jk3dkB=/" ^e*Sa'j(r-a.l)dkB, 

= |BZ|^ J(r-a-l) / e*"^^' 



Using that b'^a — 2it1 gives 



so that 



jfo.iixfo.n 



alj^ _ j 1' ^1''2 — 



0, otherwise 



1 /■ 
J =— — / JkgdkB, and consequently 
BZ Jbz 

Since the Fourier transform of E is a linear functional, it can be directly gained from the discrete Fourier 
modes of the contributing Bloch fields Ekj, . The same holds also true for the computation of the total emitted 
power by the dipole source. 

We remark that solving for Ekj, is ill-posed, when kB corresponds to a Bloch mode. Strictly, these modes do 
not appear within the Brillouin zone when lossy materials are present. However, the numerical evaluation of the 
inverse Floquet transform can be heavily affected for kB close to a resonance. This can be cured by introducing 
a small artificial damping. The so smoothed integral can be efficiently computed by an adaptive integration 
technique as demonstrated in PoUok et al."^ 

Remark 

The usage of a small artificial damping resembles the limiting absorption principle as, for example, discussed for 
photonic crystals by Joly et ali^J There, the case of a fully periodic structure is regarded (threefold periodicity in 
3D) and the limiting absorption principle became necessary to single out the correct outward radiating solution. 
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Figure 4: Emission and far field patterns generated by single dipoles (A = 570nni, position A, polarizations 
X, y, z). Left: Emission of Bloch-periodic dipoles as function of scaled Bloch- phase vector (logarithmic scale). 
The total emission is obtained by integration over the Brillouin zone. Right: Far field patterns of the isolated 
dipoles (radiated intensity as function of emission angle, where = = corresponds to normal emission, on a 
logarithmic scale). 
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Position A[nm] 77rad,x Vrnd^y ??rad,, 

A 450 0.8756 

B 450 0.8473 

C 450 0.8330 

A 570 0.9513 

B 570 0.9397 

C 570 0.9149 

A 640 0.8962 

B 640 0.9586 

C 640 0.9386 

Table 1: Computed efficiencies for different wavelengths, dipole positions and polarizations. Position A is above 
the center of a cylinder). Position B is between the first and second cylinder in x— direction x' = 500 nm, y' — 
0.0 nm), and C is the center of 4 neighboring cylinders (x' = 500 nm, y' = 500 nm). 



The same is needed here, when the twofold periodic OLED structure supports an undamped Bloch mode which 
is evanescent in the vertical direction. 

To apply the Brillouin zone integration technique to an isolated dipole J(r) = p(5(r — r') we have to solve for 
a Bloch-periodic arrangement of dipoles: 

ieN2 

Surely, due to the poor regularity it is also needed to apply the subtraction approach in this case. This can be 
done by the analytic representation of the Bloch-periodic Green's tensor, see Moroa^Sland the references therein 
for the Helmholtz equation. As an alternative, we may still use the isolated Green's function as the singular 
part. Then the correction field is no longer Bloch-periodic but jumps across the periodic boundary of the unit 
cell, 

Ec(r -f a • 1) = e**'Ba-iE^(r) + e'''S«-iE(r)^ - E^(r + a • 1). 

Fortunately, this jump condition can be seamlessly incorporated in the finite element discretization together with 
additional Neumann-type boundary terms arising in the variational form. 

5. NUMERICAL EXAMPLE 

We apply the method to a test case shown in Figure [31 Single dipoles with x-, y- and z-polarizations are placed 
at three different positions in the emitter layer of a periodically structured 3D OLED setup. The periodically 
arranged scatterers consist of silver cylinders with height 50 nm and diameter 110 nm. The grid vectors are 
Oi = [500;0;0]nm and 02 = [0;500;0]nm. Figure S] (right) shows the computed far fields for A = 570 nm and 
position A (above the center of the cylinder, x' = 0.0 nm, y' — 0.0 nm). On the left hand side the computed 
Floquet transformed total emitted power PkB,tot is shown. The emitted power for the isolated dipole is computed 
by an integration over the Brillouin zone. The sharp structures are due to the presence of complex Bloch-mode 
resonances near the Brillouin zone. To resolve this correctly an adaptive integration technique was applied. No 
artifial damping was used to mollify the integral. Table [T] gives the computed efficiencies for various wavelength 
and dipole positions. The FEM discretization was chosen to guarantee a relative error of 1% in the quantity of 
interest (?7rad)- This explains the slight asymmetry in the results (differences in 77iad,x and ?/rad,a for positions A 
and C). For generating numerical results we have used the also commercially available FEM solver JCMsuite. 

In summary, this example demonstrates that FEM based methods can accurately simulate electromagnetic 
near field distributions excited by single emitters in periodically structured media. The method can also be 
applied to arbitrary (non-periodic) structures. In post-processes, highly accurate numerical results for extraction 
efficiency of radiated power, or other derived quantities of interest can be generated. In comparison to supercell 
methods where very large computational domains have to be used (typically beyond 20 x 20 unit cells) the FEM 



computation on a single unit cell allows for very compact data space requirements and short computation times for 
single computations. The numerical integration over the Brillouin-zone can be parallelized in a straight-forward 
manner. 
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